library(magrittr)

Custom R code

It should be straightforward to use DBSLMM outputs and fam files to write custom code for CV+.

Inputs needed:

alpha <- 0.1
out <- list()
res_file <- "res-continuous.rds"
if (!file.exists(res_file)){
  res <- lapply(1:25, function(x) {
    knitr::knit_child(
      'cv-plus-child.Rmd', 
      envir = environment(), 
      quiet = TRUE
    )
    get("out")
  })
  saveRDS(res, file = res_file)
} else {
  res <- readRDS(res_file)
}
coverages <- lapply(X = res, 
                    FUN = function(mytib){
  mytib %>%
    dplyr::summarise(coverage = mean(in_interval))
       }
)
coverage_tib <- do.call("rbind", coverages) %>%
  tibble::as_tibble() %>%
  dplyr::mutate(trait_num = 1:length(res))
coverage_tib %>%
  gt::gt()
coverage trait_num
0.774 1
0.802 2
0.832 3
0.826 4
0.868 5
0.826 6
0.878 7
0.834 8
0.842 9
0.834 10
0.864 11
0.886 12
0.874 13
0.886 14
0.858 15
0.846 16
0.884 17
0.878 18
0.838 19
0.888 20
0.910 21
0.838 22
0.820 23
0.838 24
0.804 25

Boxplots of interval width

for (i in 1:length(res)){
  res[[i]] <- res[[i]] %>%
    dplyr::mutate(trait_num = i)
}
dplyr::bind_rows(res) %>%
  ggplot2::ggplot() + 
  ggplot2::geom_boxplot(ggplot2::aes(x = trait_num, y = interval_width, group = trait_num))

Incorporate estimated heritability of traits

h_tib <- readr::read_csv("~/research/ukb-intervals/shell_scripts/h2.csv", col_names = FALSE) %>% dplyr::select(-26) # drop 26, which is due to presence of trailing commas
## Rows: 5 Columns: 26
## ── Column specification ───────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (25): X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, ...
## lgl  (1): X26
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pp <- h_tib %>% 
  dplyr::summarise(dplyr::across(X1:X25, mean)) %>% 
  tidyr::pivot_longer(X1:X25) %>%
  dplyr::rename(mean_herit = value) %>%
  dplyr::mutate(trait_num = 1:25) %>%
  dplyr::select(-name) %>%
  dplyr::select(trait_num, mean_herit) %>%
  dplyr::left_join(dplyr::bind_rows(res), by = "trait_num") 
pp2 <- pp %>%
  ggplot2::ggplot() + ggplot2::geom_boxplot(ggplot2::aes(x = mean_herit, y = interval_width, group = trait_num, colour = as.factor(trait_num), label = trait_num), show.legend = FALSE)
## Warning: Ignoring unknown aesthetics: label
  pp2 %>% plotly::ggplotly(tooltip = "label")

Plot coverage vs (mean) heritability

pp2 <- pp %>%
  dplyr::group_by(trait_num) %>%
  dplyr::summarise(mean_herit = mean(mean_herit), coverage = mean(in_interval)) %>%
  ggplot2::ggplot() + ggplot2::geom_point(ggplot2::aes(x = mean_herit, y = coverage, colour = as.factor(trait_num))) 
pp2 %>%
  plotly::ggplotly()

Check correlation between (mean) heritability and coverage

summ <- pp %>%
  dplyr::group_by(trait_num) %>%
  dplyr::summarise(mean_herit = mean(mean_herit), coverage = mean(in_interval))
cor(summ$mean_herit, summ$coverage)
## [1] -0.6851291

References